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Abstract. We present a technique based on the Expectation-Maximization (EM) algorithm for the 
separation of the components of noisy mixtures in the Fourier plane. We perform a semi -blind joint 
estimation of components, mixing coefficients and noise rms levels. A priori information for the 
Qs, ■ spatial spectrum of the components and for the mixing coefficients can be naturally included in the 

f^i ' algorithm. This method is applied to the separation of distinct astrophysical emissions on simula- 

tions of future observations with the High Frequency Instrument of the Planck space mission, due to 
■ be launched in 2007. The simulations include a mixture of astrophysical emissions and instrumental 

\ white noise at the levels expected for this instrument. We have obtained good preliminary results 

O ■ with this technique, being able to blindly separate noisy mixtures with 3 components. 
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to : INTRODUCTION 

The restitution of signals or images from the observation of their mixtures has grown into 
a field of itself now classically called "source separation". Astrophysics, being a field 
of physics in which nearly all the information we can get about the physical processes 
occurring in very distant places is through observation of their electromagnetic emission, 
is naturally a field in which source separation methods can be usefully applied. One 
such application, of particular importance, can be found in millimeter and submillimeter 
astronomy. 

Mapping and interpreting sky emissions in the millimeter and sub-millimeter range 
recently made possible thanks to dedicated sensitive balloon borne and space borne in- 
struments, is indeed one of the main objectives of present and upcoming observational 
effort in astronomy. Among the scientific objectives of these observations, the precise 
measurement of primordial temperature and/or polarisation fluctuations of the Cosmic 
Microwave Background (CMB) radiation is one of the priorities, which has been given 
recently a tremendous emphasis. This radiation, emitted some 12-15 billion years ago, 
conveys a large amount of information about our universe as a whole. The importance 
of measuring anisotropics of the Cosmic Microwave Background (CMB) to constrain 
cosmological models is now well established. In the past ten years, tremendous theoreti- 
cal activity demonstrated that measuring the properties of these temperature anisotropics 
will constrain drastically the cosmological parameters describing the matter content, the 
geometry, and the evolution of our Universe [12, 13]. Recently, balloon-borne experi- 



ments such as Boomerang [4] and MAXIMA [9] have measured the CMB anisotropics in 
small patches of the sky at higher angular resolution (~ 0.2°) placing strong constraints 
on the quasi-flatness of the Universe. A new generation of satellite experiments will pro- 
vide shortly multi-frequency observations of the microwave and far infrared emission of 
the sky, with as a main objective the precise mapping of CMB fluctuations over the sky 
at high angular resolution and with unprecedented accuracy. One of these missions, the 
Microwave Anisotropy Probe, has been launched by NASA end of june 2001, and will 
provide full sky 15-30 arcminute resolution maps of the sky in three frequency channels 
with high signal to noise ratio in each pixel. Even more sensitive by an order of mag- 
nitude, the Planck mission, to be launched by ESA in 2007, will provide full sky maps 
with 5-30 arcminute resolution in 9 frequency channels between 30 and 850 GHz. 

The accuracies required for precision tests of the cosmological models, however, is 
such that it is necessary to achieve precisions on the CMB maps well below the expected 
level of contamination from astrophysical "foregrounds". Indeed, there are at least 6 
different physical emission processes which will contribute significant components in 
the Planck observations. Thus, it is crucial for the success of these future missions to 
separate CMB and foregrounds in the observed microwave maps. The separation of these 
emissions by adapted source separation methods is expected to be one of the main steps 
in the analysis of future CMB data. 

So far, two sets of independent algorithms have been proposed: MEM and Wiener 
filtering [3, 17, 11] for which the electromagnetic spectrum of the components is 
assumed known, and blind Independent Components Analysis (ICA) [1] for which no a 
priori is assumed. The former algorithms give promising results although are strongly 
limited by the uncertainties in the electromagnetic spectrum of the components which, 
as we indicated above, can be severe for some of them. The ICA algorithm has shown 
promising results for simplified non noisy mixtures but has not yet attained a sufficient 
grade of sophistication to account for instrumental noise and beam smoothing. 

We propose an alternative method for the separation of components in multi- 
frequency CMB data, based on the exploitation of the spectral diversity of the data. The 
maximisation of the likelihood is achieved with an Expectation-Maximization (EM) 
algorithm. Our method permits the simultaneous estimation of the spatial distribution 
of the components and of their electromagnetic continuum spectrum of emission. In 
section 2 we describe the basic model for noisy mixtures in the framework of the sepa- 
ration of CMB and foregrounds. In section 3 we present simulations of the HFI Planck 
observations which are used to test the separation algorithm. Section 4 describes in 
detail the EM algorithm applied to the separation of components. Section 5 summarizes 
the main results we obtain by applying the EM algorithm to our simulations. 



MODELING CMB DATA AND FOREGROUNDS 

We classify the main relevant astrophysical components in the millimetre range in three 
kinds of components. The CMB anisotropy signal, cosmological in origin, has been 
emitted in the very distant past as a relic radiation from times when the universe was 



fully ionised and before astrophysical objects as galaxies and clusters formed. Extra- 
galactic foregrounds, less distant in origin, are due to emissions coming from outside 
our galaxy. Galactic components, finally, originate from our own galaxy, and are strongly 
peaked towards the galactic plane. The main emissions at millimeter wavelengths can 
be summarised as: 

1 . CMB anisotropics. 

2. Extra- Galactic Foregrounds 

• Point sources (radio-galaxies, infrared galaxies, quasars). 

• The Sunyaev-Zeldovich (SZ) emission in clusters of galaxies. 

3. Galactic Foregrounds 

• Dust emission: thermal emission from intragalactic cold dust grains. 

• Synchrotron emission: radiation from relativistic electrons in Galactic mag- 
netic fields. 

• Free-Free (Bremsstrahlung) emission: radiation from free Galactic electrons. 

These components are known to have different spectral emission laws as a function 
of the observing frequency v. Therefore, the separation of the various emissions can be 
achieved using multi-frequency observations, i.e. the observation of the sky at different 
wavelengths, with component separation techniques based on the diversity (and possibly 
the prior knowledge) of electromagnetic spectra of foregrounds and CMB, and also on 
the spatial statistical independence of the different components. For the CMB and SZ 
effect the electromagnetic spectrum is accurately known and can be included in the 
separation methods [11]. However, for the rest of the components we dispose, in the 
best of the cases, only of spectra extrapolated from distant frequencies [5]. The spatial 
spectrum of the components is not known although reasonably good extrapolations can 
be obtained from observed data at lower resolution [3]. 

As a first step in discussing component separation techniques, we present the basic 
model which can be used to describe the observed sky emission y v [r) at position r and 
at frequency v. In the millimeter and centimeter range of the electromagnetic spectrum, 
y v (r) can be considered as a linear superposition of CMB radiation scmb{v,t) and 
foreground emissions §f(u,r) convolved with the instrumental response of the detector, 
b v (r), which is here assumed to be symmetric for simplicity. We have : 



Nf represents the number of independent foreground components considered, * defines 
the convolution operator and n v (r) is the instrumental noise of the detector. 

In the case of the CMB radiation the spatial and electromagnetic frequency depen- 
dence can be separated, 




(1) 



scmb(v,v) = gcMB{v) x s C mb{t) 






FIGURE 1. Spatial template of the CMB, dust and SZ components used in the simulations presented in 
the text. For visibility the SZ template is displayed in log scale. 



1.0000 






0.1000 






0.0100 


: CMB ■ 

rThermal SZ 




0.0010 






0.0001 







100 1000 
Frequency (GHz) 



FIGURE 2. Electromagnetic spectrum of the CMB, dust and SZ components used in the simulations 
presented in the text. These spectra fully define the missing matrix Ab when the beam smoothing effects 
are neglected. 



being gcMsiy) an d scmb(t) the CMB electromagnetic spectrum and spatial distribu- 
tion respectively. For most foreground components we can also, at least to first approxi- 
mation, separate the spatial and frequency terms [3, 11] so that equation 1 reads 

N f 

Vu(r) = 9cmb{v) s C M B (r) *b u (r) + ^g f (v) s f (r)*b u (r)+n u (r) 

where gj represents the mean electromagnetic spectrum of the foreground components. 
Typical spatial templates for CMB, dust and SZ emissions are shown in figure 1, and 
the spectral emission law (as a function of electromagnetic frequency) in figure 2. 



From the modeling point of view, the CMB component is not special and therefore 
we can write 

Vv(y,r) = ^2 9i (u) Si(r)*b v (r)+n v (r) 
i=i 

where N c = Nf + I. Thus, for a set of detectors d=l,..,N d operating at electromagnetic 
frequencies v d 

Vd(r) = ^2A di Si(r)*b d (r)+n d (r) (2) 
i=i 

where A di = gi(u d ) is a N d x N c matrix called the mixing matrix. 

If we assume the instrumental response of the detectors to be spatially invariant it is 
more convenient to work in Fourier space, since the convolution in equation 2 becomes 
a simple multiplication and we obtain 

N c 

y d (k) = J2 A ~d i (k)§ i (k)+n d (k) (3) 
i=i 

where A di {k) = A di b d {k). Further, equation 3 is verified for each Fourier mode k 
independently and if the components are assumed to be stationary random variables 
then the separation can be performed also independently at each Fourier mode. 

Hereafter, to simplify the mathematical formalism in the following sections, we will 
not include the beam smoothing in the analysis, so that equation 3 in matrix notation and 
for each mode k becomes 

y(k) = As(k)+n(k) (4) 

where y, s and n are column vectors containing N d , N c and N d complex components 
respectively, and the matrix A has dimensions N d x N c . 

Note that the separation of components requires the joint estimation of the sources 
s and the mixing matrix A (or at least some of its elements). In practice, the noise 
covariance matrix, R e — E (nn T ) is not very well known and in most cases has also to 
be estimated from the observed data, y. 



SIMULATED OBSERVATIONS FOR THE PLANCK HIGH 
FREQUENCY INSTRUMENT 

The Planck satellite will measure the emission from the whole sky in the electromag- 
netic frequency range 30 to 850 GHz. The data will be taken using two independent 
instruments, the Low Frequency Instrument (LFI) and the High Frequency Instrument 
(HFI). Here, we only consider (as a starting point for devising and testing component 
separation methods) the HFI instrument which observes the sky in 6 frequency channels 



between 100 and 850 GHz. In the frequency range of the HFI instrument, the CMB 
radiation dominates at low frequencies and the dust emission at high frequencies. In ad- 
dition, we expect significant contribution from the SZ effect in clusters of galaxies. The 
rest of the foreground components present weak emission at the HFI electromagnetic 
frequencies. 

We simulate, following the model introduced in section 2, observations from 6 detec- 
tors (N d = 6) corresponding to the 6 HFI frequency channels. To reduce the computing 
time we restrict our simulations to small patches of the sky of 300 x 300 square pixels 
of 2.5 arcmin as described in [7]. Working on small maps finds also a theoretical 
justification in the fact that the spectrum of emission of the dust (for instance) is known 
to vary slightly from a region of the sky to the other. The instrumental noise for each 
detector is considered white and Gaussian of zero mean and of rms expected for the HFI 
instrument channels. We also assume no correlation between the noise of the different 
detectors. 

The simulations include 3 components (N c = 3), CMB radiation, dust emission and 
SZ effect in clusters of galaxies. The CMB component is a Gaussian randomly generated 
realisation of CMB anisotropics obtained from a theoretical spatial power spectrum cal- 
culated with CMB FAST [16]. The electromagnetic frequency dependence of the CMB 
component is the derivative with respect to temperature of a Planck Blackbody spec- 
trum, at temperature T=2.726 K. The dust component was obtained by extrapolation of 
the IRAS satellite 100 fim map (which serves as a spatial template) to the HFI frequen- 
cies, with an electromagnetic spectrum proportional to v 2 B u (T = 17.5 K). Finally, the 
spatial template for the SZ component, due to the inverse-Compton scattering of pri- 
mordial photons by hot electrons in the ionised gaz present in clusters of galaxies, was 
obtained from a simulation of the spatial distribution of the SZ comptonization param- 
eter [6]. The emission law of the SZ emission is known theoretically, and depends only 
on the frequency of observation (in the Kompaneets approximation, which we assume 
here, see [2] for details). 

The simulated observations are obtained from the superposition of the three above spa- 
tial templates with mixing coefficients set by the electromagnetic spectra of emission 
and the frequency of observation. The observations used as inputs in the separation pro- 
cedure (with noise added at the level expected for Planck HFI detectors) are shown in 
figure 3. 



THE SEPARATION METHOD 

The blind separation of sources constitutes a common example of missing data problem 
(problems for which only partial information is available) in the literature [14]. Indeed, 
following equation 4 and considering the estimation of the mixing matrix, A, from a set 
of observations y, we can consider the components of the mixture, s as the missing data 
in our problem because they are not directly available to the observer. As discussed in 
the previous sections, we are also interested in the estimation of the noise covariance 
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FIGURE 3. Simulated observations of the six HFI channels. 



matrix, R, which can be considered as an additional parameter in the separation. 

Given the observations y, our approach is thus to estimate first the mixing ma- 
trix A and the noise covariance matrix R e by maximizing the log-likelihood function 

log p(yi.. K | A, R e ), 

(A, R e ) = argmax log p(yi.. K I A R e ) 

A. R, 

The maximization of this log-likelihood will be made with a specific implementation 
of the Expectation Maximization (EM) algorithm introduced by Dempster-Laird-Rubin 
in [8]. The estimation of the source templates is done afterwards by inverting the linear 
system y = As + n using a classical inversion method as done by Bouchet and Gispert 
in [3]. ' 



The EM algorithm: formalism 

Let's consider a missing data problem with observed data (?/(£0)fce{i..A'} an d missing 
data (a;(A;))fc e {i..A'}, for which we want to estimate the set of parameters, 0. Hereafter, 
we will note as the incomplete-data problem that for which only the observed data are 
used to obtain the unknown parameters, and the complete-data problem, that for which 
the missing data are also considered in the solution. 



An optimal estimate of 0, can be obtained by maximizing its incomplete log- 
likelihood function, 

£(0)=logp(i/|0) 

However, in many practical cases, this function is complex and difficult to work with 
and it is more convenient to solve the complete-data problem. 

The EM algorithm permits to solve the incomplete data problem. It is an iterative 
algorithm which generates a sequence of approximations to find the maximum observed 
likelihood estimator when only partial information is available, by marginalizing at each 
iteration over the missing data. 

At iteration j + 1 of the algorithm we can write 

0i+^ = M(e i ) (5) 

where M is a. mapping function, named the re-estimation transformation. After the 
initialization to an arbitrary point 0° e ©, the new estimate of is computed using 
equation (5), until a fixed point is obtained such that j+1 = j . The mapping M is 
performed in two steps 

• E-step: Computation of Q{0, j ) = E [logp(j/, x\0)\y, j ] 

• M-step: Find j+1 = argmax Q(0,0 j ) 

oe® 

where p(y\0) and p(y,x\6) denote the incomplete and the complete probability 
distributions respectively. 

A fundamental property of the EM algorithm is the fact that it ensures the 
monotonous increasing of the incomplete likelihood function. Any value of 
such that Q(9,0 j ) > Q(0 j ,0 j ) increases as well the incomplete log-likelihood, i.e, 
C(0) > Ci(0 ] ). Moreover, is a critical point of the incomplete likelihood p{y\0) 
if and only if it is a fixed point of the re-estimation transformation, M.. A more de- 
tailed description of the convergence properties of the EM algorithm can be found in [8]. 

Basic steps in the EM algorithm 

The implementation of the EM algorithm starts by choosing an initial guess for the 
unknown parameters, 0°, which is used in the first iteration of the algorithm. Then, at 
each iteration j + 1 the following basic steps are performed 

. i) e = W . 

• ii) E-step: Compute Q(0',0). 

• iii) M-step: Find j+1 such that Q(0 j+1 ,0) > Q(0',0) for all 0' e 0. 

The iterative procedure is stopped when a fixed point is reached so that j+1 = j+1 . 



Penalized EM algorithm 



In a Bayesian framework, we can consider the parameters as random variables dis- 
tributed according to an a priori p(6 | £). The a posteriori distribution for 6, p(0 | y, £) 
can be considered as a penalized version of the likelihood p(y \ 0) such that: 

p(0|y,O«p(!/|0)p(0|O 

The EM algorithm can be extended to a penalized version that converges to a local 
maximum of the penalized version of the likelihood function [10]. The penalized version 
of the EM functional Q p (0, 6 k ) is given by: 



Q P (0,0 j ) = E 
= E 



\o g p(y,x,e\o\y,e^ 

\ogp(y,x\O)\y,O>]+\ogp(0\O 



= Q M L(0,0J) + \ogp(0\0 



Application of the EM-algorithm to source separation 

Missing information model 

Now we consider the application of the EM algorithm to the problem of source 
separation in Fourier space. For this reason, we call it Spectral EM algorithm. 

As discussed in section 2, at each mode k, the vector of observations y arises from 
the following model 

y(k) = As(k) + n(k) 

Matrix A is the unknown (or partially unknown) mixing matrix and we assume that 
n{k) is a zero mean Gaussian white noise with unknown diagonal covariance matrix 
R t . We will consider = (A,R e ) the unknown parameters to be estimated. 

The incomplete likelihood function we want to maximize, p(yi..K \ A, R e ), can be 
derived by marginalizing the joint distribution p(yi..K, s±..k | A, R e ) over sources as 
follows: 

p(yi..K | A, R t ) = / p(yi..K, si.. K \ A, R e )dsi„ K 

In our method, we assume a prior knowledge of the sources spatial power spectra, i.e. 
the components s(k) are assumed to follow an a priori probability distribution of the 
form 

p(«(fc)|C(fc)) 

where (C k ) keI represent the a priori spatial power spectrum shape. 

Assuming stationarity for the sources, the Fourier modes are not correlated and we 
can perform the separation at each Fourier mode independently. The independence in 
the spectral domain can be written as: 

' p( (s(fc)W I (C(fc)W ) = \{v{s{k) I C{k)) 

kei 

< 

p ( (n(fc)W | Re ) = n p(n(fc) | 

kei 



Because of this independence, k is replaced by the integer index k and I = {1..K} is a 
fixed arbitrary arrangement, where K is the number of Fourier modes. 
We can write the complete data (yi..K, Si..k) probability distribution as 

K K 

p(yi..K, Si..k I A, R e ) = Y[p(Vk I Sk A, R e ) JJp(s fc | C k ) 

k=l k=l 

so that, the complete log-likelihood function C c is given by 
£ C (A, R e ) = \ogp(yi_,K, Si,_k | A, R € ) 

K K 



^logp(?/ fc | s k A, R e ) + ^logp(s fc | C k ) 

k=l k=l 
K 



(6) 



= 5Z(-^l 27r R A-\{Vk- As k ) 7 X 1 (Vk- As k )j + cst 

k=i v - / 

Implementation of the Spectral EM algorithm 

The functional Q(0, 9 k ) is given by 

Q(e\e j ) = E s ^[\ogp{y l .. K ,s l .. K \e)\y l .. K ,e^ 

= J si k l0gp(2/i.. K , S±.. K I 0)p(s 1 .. K I yi.. K , d^dSL.K 

and using equation 6, we get, 

Q (9 | 6 j ) = -jhg2irR € - y Tr ^(Ryy - AR T ys - R ys A T + AR SS A T )^ + est 
where 



Ryy — xY^kykUk 

Rys = -kEkVkV [s T k \y k ^] (7) 
Rss = ^Efc E [sksl\y k ,0 3 \ 

To obtain the parameter 6 j+1 = (A j+1 , R> e +1 ) at iteration j + 1, we solve the following 
gradient equations with respect to A and R e in order to maximize the functional Q. 
The complete log-likelihood function C c of equation 6 is quadratic as a function of 
A and consequently, the functional Q(A, R e \ A\ R{) is also quadratic in A and its 

derivative with respect to R e can be easily calculated. Therefore, A j+1 and R e are 
easily derived. 

§^ = KR; l (R ys -AR ss ) = 

(8) 

= f R 7 1 ( R vv ~ AR ys - R ysA T + AR SS A T )R~ 1 - fR' 1 = 



which yields to the re-estimation equations 



A j+1 — Ry S (R ss )-' 

(9) 

R{ +1 = R yy - A^R T ys - Ry S (Ai+Y + A j+1 R ss (A j+1 ) T 

We now focus on the computation of the statistics (7) involved in the re-estimation 
transformation (9). The statistic R yy , which represents the auto-co variance of the 
observations, is directly computed from the data y\..K, and remains constant throughout 
iterations of the EM algorithm. The statistic Ry S represents the a posteriori cross- 
covariance of the observations and the sources and the statistic R ss represents the a 
posteriori auto-covariance of the sources. The statistics R ys and R ss need, respectively, 
the computation of the a posteriori expectation E [s k \ y k , d j ] and the a posteriori 
covariance E [s k sjT | y k) 6 j ~\ which depends on the a priori distribution of the sources 
p(s(k) | C(k)). We consider in the following two different prior distributions: 

(i) Gaussian prior which corresponds to the Wiener filtering 

(ii) Entropic prior which corresponds to the MEM deconvolution. 

Gaussian prior: 

The a priori distribution of the source vector s k is assumed to be a zero mean Gaussian 
with covariance C k : 

p(s k | C k ) oc exp 

Then, the a posteriori distribution is 

p(s k \y k ,A,R e ) oc exp [-$(**)] 

with 

1 1 
$(s fc ) = -(y k -As k ) T R- 1 (y k -As k ) + -slC k - l s k 

Thus, the expectations are easily derived as follows 

E [s k \ =e k = [A T R^A + C~ l ] 1 A T R; 1 y k 

Eis k si] = [A T R^A + C k - 1 }- 1 + e k eT 
Entropic prior: 

For most foreground components the Gaussian prior does not fully represent the 
underlying physical processes. For this reason, Hobson et al. [11] introduced an entropic 



1 j^i-i D 

~2 k k k 



prior for the sources. Defining h k a hidden vector such that s k = L k h k , where L k is 
obtained from the Cholesky decomposition of the spectrum of sources C k , the entropic 
prior can be expressed as follows 

p(h k \m k , a) ocexp [aS(h k ,m k )] 

where S(h k) m k ) is the complex cross-entropy of h k and m k is a default model for h k 

S(h k ,m k ) = YTi=i ^kj - 2m kj - h kj log ' 



2m kj 



ki 



m k = 



Then, the a posteriori distribution of h k is: 

p(h k \y k ,A,R t ) ocexp 

with 

®( h k) = \(Vk- ALKkf 'R' 1 ^- ALh k ) - aS(h k ,m k ) 

The function <&(h k ) being convex, the a posteriori expectation of h k can be approx- 
imated by the minimum h k of this function and the covariance by the inverse of the 
Hessian [i3"(/i fc )] _1 computed at this minimum: 

E[s k ] =Lh k 

E [s kS T] = LiHihk)}- 1 ^ + E [s k ]E [s k ]* 

The minimum h k can be computed with an iterative minimization algorithm. Expres- 
sions for the gradient and the hessian of $ are easily derived, 



V hk $ = L T A T R-\ALh k -y k )+a log 



2muj 



j=l..n 



VI $ = L T A T R^ 1 AL + a diag [2m k j/ tp k j] 



APPLICATION TO PLANCK HFI SIMULATED OBSERVATIONS 

Prior information 

The penalized form of the spectral EM algorithm permits to naturally include a priori 
information about the parameters we want to estimate. In the case of CMB observations, 
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FIGURE 4. Recovered spatial template of the CMB, dust and SZ components. For visibility, the 
recovered SZ is displayed in log scale. 



one will naturally dispose of additional information from theoretical predictions or 
from previous observations at the same or other electromagnetic frequencies (although 
in most cases at lower resolution). We can take those into account in the separation 
procedure. In particular, the electromagnetic spectrum of the anisotropics of the CMB 
radiation is accurately known at centimetre and millimetre wavelengths, and therefore 
there is no need to estimate it. Further, the emission from dust grains has also been 
measured by the IRAS satellite at 100 fim providing a first class template of it. 

For a first test of the convergence capabilities of the algorithm we have only consid- 
ered weak a priori information about the spatial spectrum of the components. For the 
CMB component, the a priori spatial power spectrum was taken to be the theoretical 
model fitting best current observations [7] and for dust, the mean spatial spectrum mea- 
sured by the IRAS satellite at 100/mi [15]. A circular spatial spectrum was assumed 
as a prior for the SZ component. This spectrum was obtained from an empirical fit to 
numerical simulations of SZ templates [6]. No prior information was assumed for the 
electromagnetic spectrum of the components, i.e., no prior on the mixing matrix. 



Results 

The spectral EM algorithm has been applied to the simulated HFI observations 
presented in section 3. The prior information on the spatial spectra of the components 
is used for a Wiener filtering estimation of the missing data at each iteration of the 
algorithm. The algorithm has been implemented in C and needs about 20000 iterations 
to attain complete separation running on a 1GHz PENTIUM IV PC, for a total running 
time of about 8 hours of CPU time on 300 x 300 pixels maps. 

The recovered spatial distributions of the components are shown in figure 4. CMB 
and dust components are recovered satisfactorily, with signal to noise ratio of about 5, 
better than the signal to noise ratio in any of the detectors. The SZ is recovered with 



1.0000 








DUST^ 




0.1000 














: CMB 


. 4 






0.0100 


rTherm 


z\ SZ /^^^ 








<y 


° \ 




o \ 




0.001 


=o / 








o 


0.0001 













100 1000 
Frequency (GHz) 

FIGURE 5. The figure shows the recovered electromagnetic spectrum of the components (diamonds) 
overplotted to the input spectrum (solid lines). 



signal to noise of 1, not very satisfactorily (although the brightest clusters can be picked 
by eye on the output). This is essentially due to its peaky structure for which a Fourier 
treatment with very limited spectral a priori information is not very appropriate. 

The recovered electromagnetic spectra are shown in figure 5. The CMB and dust spec- 
tra are recovered to better than 1 % for all detectors where the level of the components 
is not negligible. The recovered electromagnetic spectrum of the SZ emission was mul- 
tiplied by a global calibration factor due to a mismatch between the assumed power 
spectrum and the true (empirical) one. In this way the true spectrum is recovered to 
better than 10% for those channels with significant SZ contribution. 



CONCLUSIONS 

The spectral version of the EM algorithm presented in this paper constitutes a useful 
tool for the semi-blind separation of astrophysical components in noisy mixtures and 
in particular to separate CMB and foreground components in future space observations 
of the CMB. A first test of the algorithm on simulated Planck HFI observations with 
3 components, CMB, dust and SZ emissions, has been proven successful. We have 
been able to achieve complete separation with signal to noise of 5 and also to recover 
the electromagnetic spectrum of the components to better than 1 % in about 20.000 
iterations. 

In the present analysis we do not consider the smearing out of the observations due to 
the finite resolution of the detectors. However, the spectral EM algorithm can be easily 
modified as proposed in section 2 to account for this. A modified version of algorithm 
is already under testing and will permit the separation of multi-frequency observations 
at different resolutions. 



Further development of the algorithm is in progress including better use of the avail- 
able prior information. In real experiments, one would like to take advantage of the fact 
that the electromagnetic spectrum is known for some of the components in which case 
there is no need to estimate it within the separation algorithm. The EM algorithm as 
presented here can naturally account for this fact by either reducing the number of free 
parameters to estimate or including a penalization term when only partial information is 
available. A first version of the algorithm implementing the former is already available 
and produces good preliminary results. 

The major drawback of our present implementation of the spectral EM algorithm is 
the computation time needed for convergence (about 20000 iterations and 8 hours of 
CPU to converge to complete separation). The acceleration of the algorithm is needed 
before it can be applied to full sky data sets or tested statistically by Monte-Carlo 
methods. Acceleration techniques are currently being explored. 
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